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Abstract. - We study the resistive transition in Josephson-junction arrays at f — 1/2 flux 
quantum per plaquette by dynamical simulations of the resistively-shunted-junction model. The 
current-voltage scaling and critical dynamics of the phases are found to be well described by 
the same critical temperature and static exponents as for the chiral (vortex-lattice) transition. 
Although this behavior is consistent with a single transition scenario, where phase and chiral 
variables order simultaneously, two different dynamic exponents result for phase coherence and 
chiral order. 



Phase transitions in two-dimensional Josephson junction arrays (JJA) has been a subject 
of much investigation [Tll!flEfllllHI?Jll71lglEfl lTUlira 

E>l|23l2Zll2EllSi]- Such arrays can be artificially fabricated [HHESE] and are also closely 
related to superconducting wire networks Experimentally, the resistive transition has 

been the one most extensively studied 2,3,4,5,6 , while theoretically several studies of XY 
models 00[i[Il[n][!]3[CIl[IlE3[H which describe the JJA, 

have been done. A significant understanding of these systems has already been achieved 
by comparing the experiments with the theoretically predicted equilibrium critical behavior 
with and without an applied magnetic field. However, to a large extent, dynamical critical 
behavior remains much less explored, particularly in the frustrated case (finite magnetic field). 
It is well known that while static critical phenomena depend on the spatial dimensionality as 
well as on the symmetry of the order parameter, the dynamic universality class will depend 
upon additional properties which do not affect the statics as, for example, conservation laws 
for the order parameter |3U) . Thus, testing the universality hypothesis of dynamical critical 
behavior requires the study of specific dynamical models. In JJA, the physically relevant 
dynamical model for the phase dynamics has not been unambiguously identified 26 , 271 1161 
117) . One would expect that, at least for an array of ideal tunnel junctions, the Resistively- 
Shunted- Junction (RSJ) model of current flow between superconducting grains would be a 
more physical representation of the system |2*8I. In experiments, the resistive transition in 
JJA is usually identified from the behavior of the current-voltage (I-V) characteristics near 
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the critical temperature. The divergent correlation length determines both the linear and 
nonlinear resistivity sufficiently close to the transition. To interpret the data and determine 
the underlying equilibrium transition, the scaling theory of Fisher et al. |31| has been widely 
used. For JJA at zero magnetic field, the resistive transition is in the Kosterlitz-Thouless 
(KT) universality class fTfYf 28I2*9| . Studies of the critical dynamics, either with Monte Carlo 
(MC) dynamics |15j or with RSJ dynamics |17| . find a behavior consistent with the dynamical 
theory of the KT transition. The exponent of the current-voltage (I-V) relation, V ~ I a , at 
the transition, assuming the universal value a = z + 1 = 3, corresponds to a dynamic exponent 
z = 2 in the resistivity scaling theory |31) . 

However, in frustrated Josephson-junction arrays (FJJA), corresponding to / = 1/2 flux 
quantum per plaquette, besides the phase variables, the vortex-lattice induced by the external 
field introduces and additional Ising-like order parameter, the chirality [7], which measures the 
direction of local current circulation in the array. The ground state consists of a pinned com- 
mensurate vortex-lattice corresponding to an antiferromagnetic arrangement of chiralities and 
vortex-lattice melting corresponds to the chiral order-disorder transition. As a consequence, 
two distinct scenarios for the phase transition have been proposed |8ll?)l ll()llTTlll2lll3lll4j : sep- 
arated phase-coherence and chiral transitions or a single transition where phase and chirality 
order simultaneously in a different universality class. Since the resistive transition corre- 
sponds to the onset of phase coherence, the later scenario where the critical dynamics should 
involve coupled variables may lead to important consequences for the resistivity scaling near 
the transition that can be detected experimentally. Previous numerical studies of the I-V 
characteristics, obtained either with RSJ dynamics |19H2l)j or with MC dynamics [22; were 
performed only for small system sizes (L < 16); while other works have studied the short- 
time dynamics of chirality , and the non-equilibrium transitions at large currents |22 | . In 
particular, the studies with RSJ dynamics used free boundary conditions to impose a driving 
current. This leads to significant additional dissipation due to boundary effects specially 
in small system sizes. In this letter, we study in detail the critical dynamics and resistivity 
scaling in FJJA by numerical simulation of the RSJ dynamics with periodic boundary condi- 
tions |23U24j including large systems sizes (L = 64, 128). This allowed us to find the following 
remarkable and unexpected results for the dynamical properties: (i) the current-voltage be- 
havior is well described by a resistive transition corresponding to the chiral transition; (ii) 
two different dynamic exponents, zxy ~ 1 and z c h ~ 2, are found for phase and chiral vari- 
ables, respectively; and (iii) at the transition, the exponent of the I-V power-law, V ~ I a , is 
a = zxy + 1 ~ 2 rather than a — 3 as for the unfrustrated case. We discuss an interpretation 
of these results within the single transition scenario including a multicritical point and also 
suggest possible explanations within the double transition picture. 

We consider a square two-dimensional array described by the overdamped RSJ model 
with current conservation at each node [2E] ■ The equations of motion for the phases 9i of the 
superconducting order parameter located at node i of the square lattice can be written as 

° j j 

where R is a uniform shunt resistance, r\ij (t) is a thermal noise with correlations (r/ij (t)rjki (£')) = 
2UbT/ R 5ij t kiS(t — t'), I is the junction critical current, Aij = f? A.dl is the line integral 
of the vector potential, and = 2wf, with the sum taken around each elementary pla- 

quette of the lattice. Dimensionless quantities are used with time in units of r = H/2eR J , 
current in units of I a , voltages in units of R I and temperature in units of HI /2ekB- A 
total current / is imposed uniformly in the array using fluctuating periodic boundary condi- 
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Fig. 1 - (a) Nonlinear resistivity E/J as a function of temperature for system size L — 64. (b) 
Scaling plot of the data for the smallest current densities. Open symbols correspond to L = 64 and 
filled ones to L = 128. 



tions |23U24| with current density J = I/L, where L is the system size and the average electric 
field E is obtained from the voltage V across the system as E = V/L = (h/2e)(d&/dt), with O 
the global phase difference or twist [231 • We integrate the dynamical equations with a second 
order Runge-Kutta-Helfand-Greenside method with time step At = 0.07r, taking 10 6 steps for 
the lowest currents while for large currents 2 x 10 4 steps were enough for proper equilibration. 
The results were averaged over 5 — 10 different initial configurations of the phases and sizes 
of L = 8 to L = 128 were considered. 

Fig. 1(a) shows the temperature dependence of the nonlinear resistivity E/J near the chiral 
transition temperature, estimated previously from equilibrium Monte Carlo simulation 8 , 
T c h — 0.455. Qualitatively, the linear resistance Rl = limj^o E/J, tends to a finite value at 
high temperatures but extrapolates to very low values at lower temperatures, independent of 
system size, consistent with the existence of a resistive transition in the range T c = 0.44 to 
0.46. In the double transition scenario, where the phase-coherence transition is expected to be 
in the KT universality class, the estimate of the KT critical temperature Tkt = 0.446 from 
Monte Carlo simulations i s very close to T c h and so without further analysis the true 
resistive critical temperature, at J = current drive, could be consistent with both estimates. 
However, we note that the IV curve for T c h > T = 0.45 > Tkt tends to zero resistivity for 
J — > 0, while the IV curve for T = 0.46 > T c h tends to finite resistivity for J — > 0. This 
already suggests that the resistivity transition occurs at T c h rather than at Tkt- In any case, 
the asymptotic critical behavior can be inferred more adequately from a scaling analysis of the 
nonlinear resistivity. According to the scaling theory |31| , measurable quantities scale with the 
diverging correlation length £ and the relaxation time r cx £ z , near the transition temperature, 
where z is the dynamical critical exponent. Then, the nonlinear resistivity should satisfy the 
scaling form 

Tj=C z 9±(^) (2) 

in two dimensions, where the + and — correspond to the behavior above and below the 
transition, respectively. For a transition in the KT universality class, the correlation length 



4 



EUROPHYSICS LETTERS 



should diverge exponentially as £ oc exp(6/|T/T c — lj 1 / 2 ), while otherwise a power-law behavior 
is expected £ oc \T/T C — 1| — with an exponent v to be determined. Thus, a scaling plot 
according to Eq. {2J can be used to verify the dynamic scaling hypotheses and the assumption 
of an underlying equilibrium transition. Such scaling plot is shown in Fig. 1(b), in the 
temperature range closest to T c h and smallest current densities, assuming the correlation 
length £ has a power-law divergence with T c = T c h and using v and z as adjustable parameters 
so that the best data collapse is obtained. Similar scaling analysis assuming a KT correlation 
length and fixing T c at the estimates of T c h or Tkt do not result in a good data collapse 
|33| . From this scaling analysis, we estimate v = 0.9(1) and the dynamical critical exponent 
z = 0.9(2). We note that the static exponent v is consistent with estimates of the chiral 
transition from equilibrium Monte Carlo simulations [5]. However, in simulations with MC 
dynamics |21| a dynamic exponent z ~ 2 was found, but this may correspond to a different 
dynamics and also note that small systems with L — 8 — 14 were analyzed. As it is well known, 
finite-size effects are very important sufficiently close to the transition when the correlation 
length £ reaches the system size L. Previous simulations with RSJ dynamics |19l20j considered 
very small sizes (L = 6 — 16). In our case, as shown in the Fig. 1(b), the two largest system 
sizes L — 64 and L — 128 give the same data collapse and so finite size effects are not dominant 
for this range of temperatures and current densities in our data. Although the above scaling 
analysis for large system sizes already suggests that the resistive transition temperature T c is 
very close to T c h with z < 2, in absence of a completely satisfactorily determination of T c from 
static critical behavior |1(JII11||T5| . from now on, we will assume j3J T c = T c h and explore to 
which extent this give us consistent results from dynamics, including finite-size effects. At T c , 
the correlation length £ will be cut off by the system size in any finite system. From Eq. , 
the nonlinear resistivity at T c should then satisfy the scaling form 

Tj = L-'g(^L) (3) 

We have tested the scaling form of Eq.© at T c = T c h for different sizes L = 8 — 128 and we 
found a very good finite size scaling with the same dynamic exponent z — 0.9(1) as shown in 
Fig. 2(a). On the contrary, when the FV curves were calculated at T = Tkt EH], it was not 
possible to obtain a reasonable data collapse using Eq.J^J for z values in the range [0.5,4]. 

It could still be argued that our estimate of z is based on a scaling analysis of the nonlinear 
I-V characteristics, which is a nonequilibrium property, and so does not reflect the underlying 
equilibrium transition. However, additional equilibrium calculations of the linear resistance 
at give a consistent estimate. From Eq. lj3|l. the linear resistance Rl = Wvcij^qE / J at 
T c should scale as Rl oc L~ z . The linear resistance can be obtained from the Kubo formula 
of equilibrium voltage fluctuations as 

Rl = (1/2T) J dt(V(t)V(0)), (4) 

without an imposing driving current, and can also be conveniently determined from the long- 
time fluctuations of the phase difference across the system |25|. Fig. 2(b), shows the finite size 
behavior of Rl- A power-law fit gives z = 0.89(6) which is in fact consistent with the estimate 
from the I-V scaling and suggests therefore that the value of z corresponds to the underlying 
equilibrium dynamical behavior. For comparison, it is also shown in Fig. 2(b) the behavior 
for the unfrustrated case, / = 0. In this case the resistive transition is in KT universality class 
and a dynamical exponent z = 2 is expected, independent of the particular dynamics [13nT7| . 
Indeed, for / = 0, the same power-law fit at the critical temperature T c — 0.887 estimated 
previously from Monte Carlo simulations |3SI gives z = 2.0(1). 
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Fig. 2 - (a) Scaling plot of the nonlinear resistivity E/J at T c — T c h = 0.455 for different system 
sizes L. (b) Linear resistance as a function of system size at the critical temperatures T c — T c h for 
/ = 1/2 and T c = 0.887 for f — 0. Power-law fits give estimates exponent z. 



Finally, to further verify that the estimate of z does in fact reflect critical phase fluctuations 
near the transition, we have also performed equilibrium calculations of the phase autocorrela- 
tion function Cxy{t) = (S(0)-S(t)), where S = Si and s = (cos(0), sin(0)). The relaxation 
time txy can be obtained from the exponential decay Cxy{t) oc exp(— t/rxy) at long times. 
Similar calculations were also performed for the chirality autocorrelation function C c h(t) to 
obtain r c fc, with the local chirality denned as \ = 2<y>(^i — @j ~ Aij)/2ir, where the sum- 
mation is taken over the elementary plaquette of the lattice and the gauge-invariant phase 
difference is restricted to the interval [— ir, w]. From dynamic finite-size scaling, the relaxation 
time should scale at T c as r oc L z , from which the z can be estimated from the slope in a 
loglog plot. Fig. 3 shows the finite-size behavior of the relaxation time at T c h for the phases 
and chiralitics. From a power-law fit we obtain zxy — 1-1(1) from the phase relaxation time 
txy which is indeed consistent with the estimate of z from the resistivity scaling discussed 
above. Naively, if the two transitions happen at the same temperature, one would expect 
that the same dynamic exponent should also hold for the chiral relaxation time. Surprisingly, 
however, the estimate from the chiral relaxation time in Fig. 3 is quite different, z c h — 2.1(1). 
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Fig. 3 - Finite-size behavior of the phase and chiral relaxation times, txy and r c ^ respectively, at 
the critical temperature T c — T c h . Power-law fits give estimates of the dynamical exponents zxy and 
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Nevertheless, this value is consistent with the result of short-time dynamics obtained from 
MC simulations [T%| . 

Thus, our numerical results imply a single transition at / = 1/2 with two dynamic expo- 
nents, zxy ~ 1 and z c h ~ 2. It should be noted that two dynamic exponents at the transition 
does not necessarily imply a breakdown of dynamic scaling in the restricted sense since 
only the chirality is an order parameter which develops true long-range order below the transi- 
tion. In two-dimensions, the phases can only develop quasi-long-range order. Since z c h > zxy, 
the longest relaxation time is still determined by the order parameter. However, extended dy- 
namic scaling 30 still applies to the resistivity and a different dynamic exponent is possible. 
Whether this only holds for the dynamics of RSJ model studied here or also for other dy- 
namical models will require further work. However, we should mention that different dynamic 
exponents for coupled order parameters have already been found previously at multicritical 
points in magnetic systems |37) . This would suggest that a possible explanation for two dy- 
namic exponents at the transition of the FJJA may rely on the existence of a multicritical 
point in the phase diagram of the relevant effective model describing the transition. Interest- 
ingly enough, the coupled XY-Ising model |Hj which is expected to describe the static critical 
behavior of the FJJA, does have a multicritical point and could be a useful framework for 
future investigations of the dynamical universality class of FJJA. 

Although the single transition scenario provides a consistent interpretation of our data, we 
should emphasize that the alternative separated transitions scenario can not be ruled out. In 
fact, recent work by Korshunov |14| argues for this scenario on good theoretical grounds. If this 
picture turns out to be the correct one then, we believe, there are two possible explanations 
for our findings: 1) the KT transition is actually much closer to T ch than estimated previously 
and so the transitions can not be resolved within the accuracy of our data; 2) the scaling 
theory of is not valid for the present case and should therefore be enlarged to include 
the interplay of two divergent length scales at nearby temperatures (131114) which can lead 
to crossover effects at small length and time scales. However, it remains unclear how the 
resistivity scaling can be so well described by the chiral transition temperature including data 
for different temperature and systems sizes. 

In conclusion, we find that the resistivity scaling and critical dynamics of FJJA are well 
described by the critical temperature corresponding to the chiral (vortex-lattice) transition. 
Two dynamic exponents, zxy ~ 1 and z c h ~ 2, are found for phase-coherence and chiral 
order, respectively, and, at the transition, the exponent of the I-V power-law, V ~ I a , is 
a = zxy + 1 ~ 2 rather than a = 3 as for the unfrustrated case. One implication of these 
results for transport experiments is that the usual method of locating the critical temperature 
from the value corresponding to a nonlinear I-V exponent a = 3, will lead to a significant 
underestimate 36 . It is worth mentioning that some experiments in overdamped arrays 3^ 
find a value of a w 2 at T c in agreement with our results. In addition, resistivity scaling 
of experimental data on wire networks [H] is also consistent with power-law correlation as 
found here but with z ~ 2. Since the dynamics of wire networks are different than the RSJ 
dynamics, this value of z may be the result of a different dynamical universality class. Further 
detailed I-V measurements combined with magnetic properties, which could in principle probe 
the chiral transition, are needed to test our results. 

* * * 
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